suppressWarnings(library(R2HTML))
suppressWarnings(library(lattice))
suppressWarnings(library(mvtnorm))

# retrieve args
Args <- commandArgs(TRUE)

#Read in data
statdata <- read.csv(Args[3], header=TRUE, sep=",")

#Copy Args
DoseResponseType <- Args[4]
ResponseVar <- Args[5]
ResponseTransform <- Args[6]
DoseVar <- Args[7]
Offsetz <- Args[8]
DoseTransform <- Args[9]
QCResponse <- Args[10]
QCDose <- Args[11]
Samples <- Args[12]
MinCoeff <- Args[13]
MaxCoeff <- Args[14]
SlopeCoeff <- Args[15]
ECIDCoeff <- Args[16]
MinStartValue <- Args[17]
MaxStartValue <- Args[18]
SlopeStartValue <- Args[19]
ECIDStartValue <- Args[20]
Equation <- Args[21]
StartValues <- Args[22]
EquationResponse <- Args[23]
EquationDose <- Args[24]


#Removing illegal charaters


if (DoseResponseType == "FourParameter")
{
YAxisTitle<-ResponseVar
XAxisTitle<-DoseVar
}

if (DoseResponseType == "FromEquation")
{
YAxisTitle<-EquationResponse
XAxisTitle<-EquationDose
}


for (i in 1:10)
{

# Additional characters included Aug 2010 (STB)
YAxisTitle<-sub("ivs_tilde_ivs"	,"~", YAxisTitle) 
YAxisTitle<-sub("ivs_star_ivs"	,"*", YAxisTitle) 
YAxisTitle<-sub("ivs_plus_ivs"	,"+", YAxisTitle) 

YAxisTitle<-sub("ivs_sp_ivs"	," ", YAxisTitle) 
YAxisTitle<-sub("ivs_ob_ivs"	,"(", YAxisTitle) 
YAxisTitle<-sub("ivs_cb_ivs"	,")", YAxisTitle) 
YAxisTitle<-sub("ivs_div_ivs"	,"/", YAxisTitle) 
YAxisTitle<-sub("ivs_pc_ivs"	,"%", YAxisTitle) 
YAxisTitle<-sub("ivs_hash_ivs"	,"#", YAxisTitle) 
YAxisTitle<-sub("ivs_pt_ivs"	,".", YAxisTitle) 
YAxisTitle<-sub("ivs_hyphen_ivs","-", YAxisTitle) 
YAxisTitle<-sub("ivs_at_ivs"	,"@", YAxisTitle) 
YAxisTitle<-sub("ivs_colon_ivs"	,":", YAxisTitle) 
YAxisTitle<-sub("ivs_exclam_ivs","!", YAxisTitle) 
YAxisTitle<-sub("ivs_quote_ivs"	,"`", YAxisTitle) 
YAxisTitle<-sub("ivs_pound_ivs"	,"£", YAxisTitle) 
YAxisTitle<-sub("ivs_dollar_ivs","$", YAxisTitle) 
YAxisTitle<-sub("ivs_hat_ivs"	,"^", YAxisTitle) 
YAxisTitle<-sub("ivs_amper_ivs"	,"&", YAxisTitle) 
YAxisTitle<-sub("ivs_obrace_ivs","{", YAxisTitle) 
YAxisTitle<-sub("ivs_cbrace_ivs","}", YAxisTitle) 
YAxisTitle<-sub("ivs_semi_ivs"	,";", YAxisTitle) 
YAxisTitle<-sub("ivs_pipe_ivs"	,"|", YAxisTitle) 
YAxisTitle<-sub("ivs_slash_ivs"	,"\\", YAxisTitle) 
YAxisTitle<-sub("ivs_osb_ivs"	,"[", YAxisTitle) 
YAxisTitle<-sub("ivs_csb_ivs"	,"]", YAxisTitle) 
YAxisTitle<-sub("ivs_eq_ivs"	,"=", YAxisTitle) 
YAxisTitle<-sub("ivs_lt_ivs"	,"<", YAxisTitle) 
YAxisTitle<-sub("ivs_gt_ivs"	,">", YAxisTitle) 
YAxisTitle<-sub("ivs_dblquote_ivs"	,"\"", YAxisTitle) 


# Additional characters included Aug 2010 (STB)
XAxisTitle<-sub("ivs_tilde_ivs"	,"~", XAxisTitle) 
XAxisTitle<-sub("ivs_star_ivs"	,"*", XAxisTitle) 
XAxisTitle<-sub("ivs_plus_ivs"	,"+", XAxisTitle) 

XAxisTitle<-sub("ivs_sp_ivs"	," ", XAxisTitle) 
XAxisTitle<-sub("ivs_ob_ivs"	,"(", XAxisTitle) 
XAxisTitle<-sub("ivs_cb_ivs"	,")", XAxisTitle) 
XAxisTitle<-sub("ivs_div_ivs"	,"/", XAxisTitle) 
XAxisTitle<-sub("ivs_pc_ivs"	,"%", XAxisTitle) 
XAxisTitle<-sub("ivs_hash_ivs"	,"#", XAxisTitle) 
XAxisTitle<-sub("ivs_pt_ivs"	,".", XAxisTitle) 
XAxisTitle<-sub("ivs_hyphen_ivs","-", XAxisTitle) 
XAxisTitle<-sub("ivs_at_ivs"	,"@", XAxisTitle) 
XAxisTitle<-sub("ivs_colon_ivs"	,":", XAxisTitle) 
XAxisTitle<-sub("ivs_exclam_ivs","!", XAxisTitle) 
XAxisTitle<-sub("ivs_quote_ivs"	,"`", XAxisTitle) 
XAxisTitle<-sub("ivs_pound_ivs"	,"£", XAxisTitle) 
XAxisTitle<-sub("ivs_dollar_ivs","$", XAxisTitle) 
XAxisTitle<-sub("ivs_hat_ivs"	,"^", XAxisTitle) 
XAxisTitle<-sub("ivs_amper_ivs"	,"&", XAxisTitle) 
XAxisTitle<-sub("ivs_obrace_ivs","{", XAxisTitle) 
XAxisTitle<-sub("ivs_cbrace_ivs","}", XAxisTitle) 
XAxisTitle<-sub("ivs_semi_ivs"	,";", XAxisTitle) 
XAxisTitle<-sub("ivs_pipe_ivs"	,"|", XAxisTitle) 
XAxisTitle<-sub("ivs_slash_ivs"	,"\\", XAxisTitle) 
XAxisTitle<-sub("ivs_osb_ivs"	,"[", XAxisTitle) 
XAxisTitle<-sub("ivs_csb_ivs"	,"]", XAxisTitle) 
XAxisTitle<-sub("ivs_eq_ivs"	,"=", XAxisTitle) 
XAxisTitle<-sub("ivs_lt_ivs"	,"<", XAxisTitle) 
XAxisTitle<-sub("ivs_gt_ivs"	,">", XAxisTitle) 
XAxisTitle<-sub("ivs_dblquote_ivs"	,"\"", XAxisTitle) 
}






if (DoseResponseType == "FromEquation")
{
	#Setup the html file and associated css file
	htmlFile <- sub(".csv", ".html", Args[3]); #determine the file name of the html file
	.HTML.file = htmlFile
	cssFile <- "r2html.css"
	cssFile <- paste("'",cssFile,"'", sep="") #need to enclose in quotes when path has spaces in it
	HTMLCSS(CSSfile = cssFile)
	
	#Output HTML header
	HTML.title("<bf>SilveR Dose-Response Analysis", HR=1, align="left")


	#Titles and description

	title<-c("Variable selection")
	HTML.title(title, HR=2, align="left")

	add<-paste(c("The  "), EquationResponse, sep="")
	add<-paste(add, " response is currently being analysed using the user defined equation in the Dose-Response Analysis module. ", sep="")
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(add, HR=0, align="left")

	#Dose fit equation

	statdata$x = eval(parse(text = paste("statdata$", EquationDose)))
	statdata$y = eval(parse(text = paste("statdata$", EquationResponse )))
	Equation2 <- eval(parse(text = paste("y~", Equation)))



	#Seperate out the start values

	tempChanges <-strsplit(StartValues, ",")

	txtexpectedChanges <- c("")
	for(i in 1:length(tempChanges[[1]])) 
	{ 
		txtexpectedChanges [length(txtexpectedChanges)+1]=(tempChanges[[1]][i]) 
	}
	tabs<-matrix(nrow=(length(txtexpectedChanges)-1), ncol=2)

	for (i in 2:length(txtexpectedChanges))
	{
		Changes <-strsplit(txtexpectedChanges[i], "=")
		txtChanges <- c("")
		for(j in 1:length(Changes[[1]])) 
		{ 
			txtChanges [length(txtChanges)+1]=(Changes[[1]][j]) 
		}
	tabs[(i-1),1]=txtChanges[2]
	tabs[(i-1),2]=as.numeric(txtChanges[3])*1
	}

	#Setting up the list of start values

	nameparas<-c()
	paras<-c()
	index<-1
	for (i in 1:(length(txtexpectedChanges))-1)
	{
	paras[i]=as.numeric(tabs[i,2])
	nameparas[i]=tabs[i,1]
	}
	names(paras)<-nameparas

	dosefit<-nls(Equation2, start=paras, data=statdata)

if (min(statdata$x) > 0)
{
xmin<-min(statdata$x)*0.9
} else 
{
xmin<-min(statdata$x)*1.1
}

if (max(statdata$x) > 0)
{
xmax<-max(statdata$x)*1.1
} else 
{
xmax<-max(statdata$x)*0.9
}


if (min(statdata$y) > 0)
{
ymino<-min(statdata$y)*0.9
} else 
{
ymino<-min(statdata$y)*1.1
}

if (max(statdata$y) > 0)
{
ymaxo<-max(statdata$y)*1.1
} else 
{
ymaxo<-max(statdata$y)*0.9
}



av<-seq(xmin,xmax,0.01)
bv<-predict(dosefit,list(x =av))


if (min(bv) > 0)
{
yminp<-min(bv)*0.9
} else 
{
yminp<-min(bv)*1.1
}

if (max(bv) > 0)
{
ymaxp<-max(bv)*1.1
} else 
{
ymaxp<-max(bv)*0.9
}



if (ymino < yminp)
{
ymin <- ymino
} else 
{
ymin <-yminp
}

if (ymaxo > ymaxp)
{
ymax <- ymaxo
} else 
{
ymax <-ymaxp
}


	#Table of parameter estimates

	table<-summary(dosefit)$parameters

	tablen<-length(unique(rownames(table)))
	tabz<-matrix(nrow=tablen, ncol=4)
	
	for (i in 1:tablen)
	{
		tabz[i,1]=format(round(table[i,1], 2), nsmall=2, scientific=FALSE)
	}

	for (i in 1:tablen)
	{
		tabz[i,2]=format(round(table[i,2], 3), nsmall=3, scientific=FALSE)
	}

	for (i in 1:tablen)
	{
		tabz[i,3]=format(round(table[i,3], 2), nsmall=2, scientific=FALSE)
	}

	for (i in 1:tablen)
	{
		tabz[i,4]=format(round(table[i,4], 4), nsmall=4, scientific=FALSE)
	}

	for (i in 1:tablen) 
	{
		if (tabz[i,4]<0.001) 
		{
			tabz[i,4]<-0.001
			tabz[i,4]<- paste("<",tabz[i,4])
		}
	}

	header<-c(" ", " "," ", " ")
	tabs<-rbind(header, tabz)
	rownames(tabs)<-c("Parameter", rownames(table))
	colnames(tabs)<-c("Estimate", "Std Error", "t-value", "p-value")



	#DF test
	df<-length(na.omit(eval(parse(text = paste("statdata$", EquationResponse )))))-tablen
	if (df <= 4)
	{
		warning<-c("Unfortunately the residual degrees of freedom are low (less than 5). This may make the estimation of the underlying variability, and hence the results of the statistical tests, unreliable. This can be caused by attempting to model too many parameters. We recommend you fix some of the parameters.")  
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(warning, HR=0, align="left")
		quit()
	}


	#Plotting the results
	statdata$conczzzz = eval(parse(text = paste("statdata$", EquationDose)))
	statdata$respzzzz = eval(parse(text = paste("statdata$", EquationResponse )))

	HTMLbr()
	title<-c("Scatterplot of raw data including the predicted fit")
	HTML.title(title, HR=2, align="left")

	scatterPlot <- sub(".html", "scatterPlot.jpg", htmlFile)
	jpeg(scatterPlot)
	par(las=1)
	plot(eval(parse(text = paste("statdata$", EquationDose))),eval(parse(text = paste("statdata$", EquationResponse ))), xlab=XAxisTitle , ylab = YAxisTitle,  ylim = c(ymin,ymax))
	
	bv<-predict(dosefit,list(x =av))
	lines(av,bv,col="black")

	void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", scatterPlot), Align="centre")



	#Table of parameter estimates print
	HTMLbr()
	title<-c("Table of model parameters and summary statistics")
	HTML.title(title, HR=2, align="left")
	HTML(tabs, classfirstline="second", align="left")
}

	
if (DoseResponseType == "FourParameter")
{

	if(unique(eval(parse(text = paste("statdata$", DoseVar))))[1] == 0)
	{	
		# Setting up the offset parameter
		if (Offsetz  == "NULL" && DoseTransform == "Log10")
		{
			offset <-(unique(eval(parse(text = paste("statdata$", DoseVar))))[2])/10
		} else 
		if (Offsetz  == "NULL" && DoseTransform == "Loge")
		{
			offset <-(unique(eval(parse(text = paste("statdata$", DoseVar))))[2])/exp(1)
		} else
		if (Offsetz  != "NULL" ) 
		{
			offset<- as.numeric(Offsetz)
		}
	} else
	if(unique(eval(parse(text = paste("statdata$", DoseVar))))[1] != 0)
	{
		# Setting up the offset parameter
		if (Offsetz  == "NULL" && DoseTransform == "Log10")
		{
			offset <-(unique(eval(parse(text = paste("statdata$", DoseVar))))[1])/10
		} else 
		if (Offsetz  == "NULL" && DoseTransform == "Loge")
		{
			offset <-(unique(eval(parse(text = paste("statdata$", DoseVar))))[1])/exp(1)
		} else
		if (Offsetz  != "NULL" ) 
		{
			offset<- as.numeric(Offsetz)
		}
	} 


# Setting up the fixed parameter


if (MinCoeff !="NULL")
{
	MinCoeffp <- as.numeric(MinCoeff)
}
if (MaxCoeff !="NULL")
{
	MaxCoeffp <- as.numeric(MaxCoeff)
}
if (SlopeCoeff !="NULL")
{
	SlopeCoeffp <- as.numeric(SlopeCoeff)
}
if (ECIDCoeff !="NULL" && DoseTransform == "Log10")
{
	ECIDCoeffp <- log10(as.numeric(ECIDCoeff)+offset)
} else 
if (ECIDCoeff !="NULL" && DoseTransform == "Loge")
{
	ECIDCoeffp <- log(as.numeric(ECIDCoeff)+offset)
}

#Setup the html file and associated css file
htmlFile <- sub(".csv", ".html", Args[3]); #determine the file name of the html file
.HTML.file = htmlFile
cssFile <- "r2html.css"
cssFile <- paste("'",cssFile,"'", sep="") #need to enclose in quotes when path has spaces in it
HTMLCSS(CSSfile = cssFile)

#Output HTML header
HTML.title("<bf>SilveR Dose-Response Analysis", HR=1, align="left")


# Setting up the concentration parameters

statdata$responsezzzz = eval(parse(text = paste("statdata$", ResponseVar)))

if (DoseTransform == "Log10")
{
	statdata$logconczzzz = log10(eval(parse(text = paste("statdata$", DoseVar)))+offset)
	statdata$conczzzz = eval(parse(text = paste("statdata$", DoseVar)))
} else  
if (DoseTransform == "Loge")
{
	statdata$logconczzzz = log(eval(parse(text = paste("statdata$", DoseVar)))+offset)
	statdata$conczzzz = eval(parse(text = paste("statdata$", DoseVar)))
}


#Setting up graph titles

if (ResponseTransform == "None")
{
	YAxisTitle <- YAxisTitle
} else 
{
	Yadd <- paste(ResponseTransform, " (", sep="")
	Yadd <- paste(Yadd , YAxisTitle, sep="")
	Yadd <- paste(Yadd , ")", sep="")
	YAxisTitle<-Yadd
}


if (DoseTransform == "Log10")
{
	Xadd <- paste(XAxisTitle, " (on the Log10 scale)",sep="")
	XAxisTitle<-Xadd
} else 
if (DoseTransform == "Loge")
{
	Xadd <- paste(XAxisTitle, " (on the Loge scale)",sep="")
	XAxisTitle<-Xadd
}


# Setting up the start parameters - need to check user defined options

if (MinStartValue == "NULL")
{
	minp<-min(unlist(lapply(split(eval(parse(text = paste("statdata$", ResponseVar))),statdata$logconczzzz),mean)))
} else 
{
	minp<-MinStartValue
}

if (MaxStartValue == "NULL")
{
	maxp<-max(unlist(lapply(split(eval(parse(text = paste("statdata$", ResponseVar))),statdata$logconczzzz),mean)))
} else 
{
	maxp<-MaxStartValue
}

if (SlopeStartValue == "NULL")
{
	slopep<-1
	if(maxp>minp) 
	{	
		slopep=1
	} else 
	{
		slopep=-1
	}
} else 
{
	slopep<-SlopeStartValue
}


if (DoseTransform == "Log10" && ECIDStartValue != "NULL")
{
	ed50p<-log10(as.numeric(ECIDStartValue))
} else 
if (DoseTransform == "Loge" && ECIDStartValue != "NULL")
{
	ed50p<-log(as.numeric(ECIDStartValue))	
} else  
if (ECIDStartValue == "NULL")
{
	temp<-sort(unique(statdata$logconczzzz))
	ed50p<-mean(temp[-1])
}


#Titles and description

title<-c("Response and dose variables")
HTML.title(title, HR=2, align="left")


add<-paste(c("The  "), ResponseVar, sep="")
add<-paste(add, " response is currently being analysed by the Dose-Response Analysis module. ", sep="")


if(ResponseTransform != "None")
{
	add<-paste(add, "The response has been ", sep="")
	add<-paste(add, ResponseTransform, sep="")
	add<-paste(add, " transformed prior to analysis.", sep="")
}

HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")


if(DoseTransform == "Log10")
{
	add<-paste(c("The dose variable ("), DoseVar, sep="")
	add<-paste(add, ") has been Log10 transformed prior to analysis.", sep="")
} else 
if(DoseTransform == "Loge")
{
	add<-paste(c("The dose variable ("), DoseVar, sep="")
	add<-paste(add, ") has been Loge transformed prior to analysis.", sep="")
}

HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")


#Scatterplot
#HTMLbr()
#title<-c("Scatterplot of the raw data")
#HTML.title(title, HR=2, align="left")
#
#scatterPlot2 <- sub(".html", "scatterPlot2.jpg", htmlFile)
#jpeg(scatterPlot2)
#xyplot(responsezzzz~logconczzzz, statdata, col="black", xlab=XAxisTitle, ylab = YAxisTitle)#, scales=list(x=list(rot=90)) )
#void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", scatterPlot2), Align="centre")
#
#HTML.title("</bf> ", HR=2, align="left")
#HTML.title("</bf> Tip: Use this plot to identify start values, if you have convergence issues.", HR=0, align="left")



#Fitting the four paramter model

if        (MinCoeff == "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff == "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~D         + (A        -D        )/(1+10^((C         -logconczzzz)*B          )), start=list(A=maxp,B=slopep,C= ed50p,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff == "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~D         + (A        -D        )/(1+exp((C         -logconczzzz)*B          )), start=list(A=maxp,B=slopep,C= ed50p,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff == "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~D         + (MaxCoeffp-D        )/(1+10^((C         -logconczzzz)*B          )), start=list(B=slopep,C= ed50p,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff == "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~D         + (MaxCoeffp-D        )/(1+exp((C         -logconczzzz)*B          )), start=list(B=slopep,C= ed50p,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff == "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~D         + (A        -D        )/(1+10^((C         -logconczzzz)*SlopeCoeffp)), start=list(A=maxp,C= ed50p,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff == "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~D         + (A        -D        )/(1+exp((C         -logconczzzz)*SlopeCoeffp)), start=list(A=maxp,C= ed50p,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff != "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~D         + (A        -D        )/(1+10^((ECIDCoeffp-logconczzzz)*B          )), start=list(A=maxp,B=slopep,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff != "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~D         + (A        -D        )/(1+exp((ECIDCoeffp-logconczzzz)*B          )), start=list(A=maxp,B=slopep,D=minp), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff == "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (A        -MinCoeffp)/(1+10^((C         -logconczzzz)*B          )), start=list(A=maxp,B=slopep,C= ed50p), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff == "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (A        -MinCoeffp)/(1+exp((C         -logconczzzz)*B          )), start=list(A=maxp,B=slopep,C= ed50p), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff == "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~D         + (MaxCoeffp-D        )/(1+10^((C         -logconczzzz)*SlopeCoeffp)), start=list(C= ed50p,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff == "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~D         + (MaxCoeffp-D        )/(1+exp((C         -logconczzzz)*SlopeCoeffp)), start=list(C= ed50p,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff != "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~D         + (MaxCoeffp-D        )/(1+10^((ECIDCoeffp-logconczzzz)*B          )), start=list(B=slopep,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff != "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~D         + (MaxCoeffp-D        )/(1+exp((ECIDCoeffp-logconczzzz)*B          )), start=list(B=slopep,D=minp), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff == "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (MaxCoeffp-MinCoeffp)/(1+10^((C         -logconczzzz)*B          )), start=list(B=slopep,C= ed50p), data=statdata)
} else if (MinCoeff != "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff == "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (MaxCoeffp-MinCoeffp)/(1+exp((C         -logconczzzz)*B          )), start=list(B=slopep,C= ed50p), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff != "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~D         + (A        -D        )/(1+10^((ECIDCoeffp-logconczzzz)*SlopeCoeffp)), start=list(A=maxp,D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff != "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~D         + (A        -D        )/(1+exp((ECIDCoeffp-logconczzzz)*SlopeCoeffp)), start=list(A=maxp,D=minp), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff == "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (A        -MinCoeffp)/(1+10^((C         -logconczzzz)*SlopeCoeffp)), start=list(A=maxp,C= ed50p), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff == "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (A        -MinCoeffp)/(1+exp((C         -logconczzzz)*SlopeCoeffp)), start=list(A=maxp,C= ed50p), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff != "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (A        -MinCoeffp)/(1+10^((ECIDCoeffp-logconczzzz)*B          )), start=list(A=maxp,B=slopep), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff != "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (A        -MinCoeffp)/(1+exp((ECIDCoeffp-logconczzzz)*B          )), start=list(A=maxp,B=slopep), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff != "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~D         + (MaxCoeffp-D        )/(1+10^((ECIDCoeffp-logconczzzz)*SlopeCoeffp)), start=list(D=minp), data=statdata) 
} else if (MinCoeff == "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff != "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~D         + (MaxCoeffp-D        )/(1+exp((ECIDCoeffp-logconczzzz)*SlopeCoeffp)), start=list(D=minp), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff == "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (MaxCoeffp-MinCoeffp)/(1+10^((C         -logconczzzz)*SlopeCoeffp)), start=list(C= ed50p), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff == "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (MaxCoeffp-MinCoeffp)/(1+exp((C         -logconczzzz)*SlopeCoeffp)), start=list(C= ed50p), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff != "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (MaxCoeffp-MinCoeffp)/(1+10^((ECIDCoeffp-logconczzzz)*B          )), start=list(B=slopep), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff == "NULL"&& ECIDCoeff != "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (MaxCoeffp-MinCoeffp)/(1+exp((ECIDCoeffp-logconczzzz)*B          )), start=list(B=slopep), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff != "NULL"&& DoseTransform != "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (A        -MinCoeffp)/(1+10^((ECIDCoeffp-logconczzzz)*SlopeCoeffp)), start=list(A=maxp), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff == "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff != "NULL"&& DoseTransform == "Loge") {dosefit<-nls(responsezzzz~MinCoeffp + (A        -MinCoeffp)/(1+exp((ECIDCoeffp-logconczzzz)*SlopeCoeffp)), start=list(A=maxp), data=statdata) 
} else if (MinCoeff != "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff != "NULL"&& DoseTransform != "Loge") 
	{

	HTMLbr()
	title<-c("Warning")
	HTML.title(title, HR=2, align="left")
	warning<-c("You need to estimate at least one of the parameters, hence no analysis has been performed.")  
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(warning, HR=0, align="left")
	} else if (MinCoeff != "NULL"&& MaxCoeff != "NULL"&& SlopeCoeff != "NULL"&& ECIDCoeff != "NULL"&& DoseTransform == "Loge") 
		{
		HTMLbr()
		title<-c("Warning")
		HTML.title(title, HR=2, align="left")		
		warning<-c("You need to estimate at least one of the parameters, hence no analysis has been performed.")  
		HTML.title("</bf> ", HR=2, align="left")
		HTML.title(warning, HR=0, align="left")
		} 



	#DF test
	table<-summary(dosefit)$parameters
	tablen<-length(unique(rownames(table)))
	df<-length(na.omit(eval(parse(text = paste("statdata$", ResponseVar )))))-tablen
	if (df <= 4)
	{
	warning<-c("Unfortunately the residual degrees of freedom are low (less than 5). This may make the estimation of the underlying variability, and hence the results of the statistical tests, unreliable. This can be caused by attempting to model too many parameters. We recommend you fix some of the parameters.")  
	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(warning, HR=0, align="left")
	quit()
	}


#Plotting the results

HTMLbr()
title<-c("Scatterplot of raw data including the predicted fit")
HTML.title(title, HR=2, align="left")

scatterPlot <- sub(".html", "scatterPlot.jpg", htmlFile)
jpeg(scatterPlot)
par(las=1)
plot(statdata$logconczzzz,eval(parse(text = paste("statdata$", ResponseVar))), xlab=XAxisTitle, ylab = YAxisTitle  )
av<-seq(-10,10,0.01)
bv<-predict(dosefit,list(logconczzzz=av))
lines(av,bv,col="black")

void <- HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", scatterPlot), Align="centre")


#Table of parameter estimates

HTMLbr()
title<-c("Table of curve parameters and summary statistics")
HTML.title(title, HR=2, align="left")

table<-summary(dosefit)$parameters
tablen<-length(unique(rownames(table)))

paratp<-rownames(table)
paranams<-c(1:tablen)

for (i in 1:tablen)
{ 
	if (paratp[i] == "A")
	{
		paranams[i] = "Max/Min"
	} else 
	if (paratp[i] == "B")
	{
		paranams[i] = "Slope"
	} else
	if (paratp[i] == "C")
	{
		paranams[i] = "logED50"
	} else
	if (paratp[i] == "D")
	{
		paranams[i] = "Min/Max"
	} 
}


tabz<-matrix(nrow=tablen, ncol=4)

for (i in 1:tablen)
{
	tabz[i,1]=format(round(table[i,1], 2), nsmall=2, scientific=FALSE)
}

for (i in 1:tablen)
{
	tabz[i,2]=format(round(table[i,2], 3), nsmall=3, scientific=FALSE)
}

for (i in 1:tablen)
{
	tabz[i,3]=format(round(table[i,3], 2), nsmall=2, scientific=FALSE)
}

for (i in 1:tablen)
{
	tabz[i,4]=format(round(table[i,4], 4), nsmall=4, scientific=FALSE)
}

for (i in 1:tablen) 
{
	if (tabz[i,4]<0.001) 
	{
		tabz[i,4]<-0.001
		tabz[i,4]<- paste("<",tabz[i,4])
	}
}

header<-c(" ", " "," ", " ")
tabs<-rbind(header, tabz)

rownames(tabs)<-c("Parameter", paranams)
colnames(tabs)<-c("Estimate", "Std Error", "t-value", "p-value")
HTML(tabs, classfirstline="second", align="left")



ictab<-data.frame(table)

if (ECIDCoeff == "NULL")
{
	#Table of parameter estimates

	HTMLbr()
	title<-c("Back transformed ED50 estimate with 95% confidence intervals")
	HTML.title(title, HR=2, align="left")

	#No. of degrees of freedom
	df<-length(statdata$responsezzzz)-tablen

	critval95<-qmvt(0.95, df = df, tail = "both", abseps=0.000001)[1]

	for (i in 1:tablen)
	{
		if (rownames(ictab)[i] == "C")
		{
			logED50<-ictab[i,1]
			ED50ERR<-ictab[i,2]

			loglower<-logED50-as.numeric(critval95)*ED50ERR
			logupper<-logED50+as.numeric(critval95)*ED50ERR	

	
			if(DoseTransform == "Log10")
			{
				ED50<-10**logED50-offset
				lower<-10**loglower-offset
				upper<-10**logupper-offset
			}		
			if(DoseTransform == "Loge")
			{
				ED50<-exp(logED50)-offset
				lower<-exp(loglower)-offset
				upper<-exp(logupper)-offset
			}
		}
	}

	ED50table<-matrix(nrow=1, ncol=3)
	ED50table[1,1]<-format(round(ED50,2), nsmall=2, scientific=FALSE)
	ED50table[1,2]<-format(round(lower,2), nsmall=2, scientific=FALSE)
	ED50table[1,3]<-format(round(upper,2), nsmall=2, scientific=FALSE)

	ED50table<-data.frame(ED50table)
	rownames(ED50table)<-c("Estimate")
	colnames(ED50table)<-c("ED50", "Lower 95% CI", "Upper 95% CI")
	HTML(ED50table, classfirstline="second", align="left")
}


#Extracting parameter estimates
if (MinCoeff !="NULL")
{
D <- as.numeric(MinCoeff)
} else { 
		for (i in 1:tablen)
		{
			if (rownames(ictab)[i] == "D")
			{
			D <-ictab[i,1]
			}
		}
	}

if (MaxCoeff !="NULL")
{
A <- as.numeric(MaxCoeff)
} else { 
		for (i in 1:tablen)
		{
			if (rownames(ictab)[i] == "A")
			{
			A <-ictab[i,1]
			}
		}
	}
if (SlopeCoeff !="NULL")
{
B <- as.numeric(SlopeCoeff)
} else { 
		for (i in 1:tablen)
		{
			if (rownames(ictab)[i] == "B")
			{
			B <-ictab[i,1]
			}
		}
	}
if (ECIDCoeff !="NULL" && DoseTransform != "Loge")
{
	C <- log10(as.numeric(ECIDCoeff)+offset)
} else 
if (ECIDCoeff !="NULL" && DoseTransform == "Loge")
{
	C <- log(as.numeric(ECIDCoeff)+offset)
} else	
if (ECIDCoeff =="NULL" )
{
	for (i in 1:tablen)
	{
		if (rownames(ictab)[i] == "C")
		{
			C <-ictab[i,1]
		}
	}
}



if (ECIDCoeff !="NULL" && DoseTransform == "Log10")
{
	btC <- as.numeric(ECIDCoeff)
} else 
if (ECIDCoeff !="NULL" && DoseTransform == "Loge")
{
	btC <- as.numeric(ECIDCoeff)
} else	
if (ECIDCoeff =="NULL" && DoseTransform == "Log10")
{
	for (i in 1:tablen)
	{
		if (rownames(ictab)[i] == "C")
		{
			btC <-10**(ictab[i,1])-offset
		}
	}
} else 
if (ECIDCoeff =="NULL" && DoseTransform == "Loge")
{
	for (i in 1:tablen)
	{
		if (rownames(ictab)[i] == "C")
		{
			btC <-exp(ictab[i,1])-offset
		}
	}
}


#Analysis decription

HTMLbr()
title<-c("Description of the analysis results")
HTML.title(title, HR=2, align="left")

add<-c("The data was analysed using non-linear regression. ")

if (MinCoeff == "NULL")
{
	add<-paste(add, "The estimate of the minimum of the curve is ", sep="")
	if (D<A){
		add<-paste(add, format(round(D,2), nsmall=2, scientific=FALSE), sep="")
		} else {
			add<-paste(add, format(round(A,2), nsmall=2, scientific=FALSE), sep="")
			}
	add<-paste(add, ". ", sep="")
}

if (MaxCoeff == "NULL")
{
	add<-paste(add, "The estimate of the maximum of the curve is ", sep="")
	if (A>D){
		add<-paste(add, format(round(A,2), nsmall=2, scientific=FALSE), sep="")
		} else {
			add<-paste(add, format(round(D,2), nsmall=2, scientific=FALSE), sep="")
			}
	add<-paste(add, ". ", sep="")
}

if (SlopeCoeff == "NULL")
{
	add<-paste(add, "The estimate of the slope coefficient of the curve is ", sep="")
	add<-paste(add, format(round(B,2), nsmall=2, scientific=FALSE), sep="")
	add<-paste(add, ". ", sep="")
}

if (ECIDCoeff == "NULL")
{
	add<-paste(add, "The estimate of the logED50 of the curve is ", sep="")
	add<-paste(add, format(round(C,2), nsmall=2, scientific=FALSE), sep="")
	add<-paste(add, ". ", sep="")
}

if (ECIDCoeff == "NULL")
{
	add<-paste(add, "The back transformed estimate of the ED50 of the curve is ", sep="")
	add<-paste(add, format(round(btC,3), nsmall=3, scientific=FALSE), sep="")
	add<-paste(add, ". ", sep="")
}

if (MinCoeff != "NULL")
{
	add<-paste(add, "The minimum of the curve was fixed at ", sep="")
	add<-paste(add, MinCoeff, sep="")
	add<-paste(add, ". ", sep="")
}

if (MaxCoeff != "NULL")
{
	add<-paste(add, "The maximum of the curve was fixed at ", sep="")
	add<-paste(add, MaxCoeff, sep="")
	add<-paste(add, ". ", sep="")
}

if (SlopeCoeff != "NULL")
{
	add<-paste(add, "The slope of the curve was fixed at ", sep="")
	add<-paste(add, SlopeCoeff, sep="")
	add<-paste(add, ". ", sep="")
}

if (ECIDCoeff != "NULL")
{
	add<-paste(add, "The ED50 of the curve was fixed at ", sep="")
	add<-paste(add, ECIDCoeff, sep="")
	add<-paste(add, ". ", sep="")
}

HTML.title("</bf> ", HR=2, align="left")
HTML.title(add, HR=0, align="left")

#QC samples

if (QCResponse != "NULL" && QCDose != "NULL")
{
	HTMLbr()
	title<-c("Scatterplot of responses and quality control (QC) samples")
	HTML.title(title, HR=2, align="left")

	# Plot of data

	statdata$QCresponsezzzz = eval(parse(text = paste("statdata$", QCResponse)))


	if (DoseTransform == "Log10")
	{	
		statdata$logQCconczzzz = log10(eval(parse(text = paste("statdata$", QCDose))))
		statdata$QCconczzzz = eval(parse(text = paste("statdata$", QCDose)))
	} else	if (DoseTransform == "Loge")
	{	
		statdata$logQCconczzzz = log(eval(parse(text = paste("statdata$", QCDose))))
		statdata$QCconczzzz = eval(parse(text = paste("statdata$", QCDose)))
	}


	# setting up the dataset

	SC<-data.frame(cbind(statdata$responsezzzz, statdata$logconczzzz,statdata$conczzzz))
	SClen<-dim(SC)[1]
	for (i in 1:SClen)
	{
		SC$Type[i]="Assay standard curve"
	}
	colnames(SC)<-c("response","logdose","dose","Type")

	QCs<-data.frame(cbind(statdata$QCresponsezzzz, statdata$logQCconczzzz, statdata$QCconczzzz))
	QCs<-na.omit(QCs) 
	QClen<-dim(QCs)[1]
	for (i in 1:QClen)
	{
		QCs$Type[i]="Quality controls"
	}

	colnames(QCs)<-c("response","logdose","dose","Type")

	alldata <-data.frame(rbind(SC,QCs))

	rows<-dim(alldata)[1]
	cols<-dim(alldata)[2]
	nlevels<-length(unique(as.factor(alldata$Type)))
	extra<-matrix(data=NA, nrow=rows, ncol=nlevels)
	for (i in 1:nlevels)
	{
		for (j in 1:rows)
		{
			if (alldata$Type[j] == unique(as.factor(alldata$Type))[i])
			{
				extra[j,i]<-alldata$response[j]
			}
		}
	}
	newdata<-cbind(alldata, extra)
	catplotdata<-data.frame(newdata)
	for (k in 1:nlevels)
	{
		tempdata<-catplotdata
		tempdata2<-subset(tempdata, tempdata$Type == unique(levels(as.factor(tempdata$Type)))[k])
	}
	index<-c(1:nlevels)
	newnames<-c(colnames(alldata),index)
	colnames(catplotdata)<-newnames
	ncscatterplot3 <- sub(".html", "ncscatterplot3.jpg", htmlFile)
	jpeg(ncscatterplot3)

	if (is.numeric(alldata$Type)==FALSE)
	{
		#Adjusting y  axis to fit in legend
		maxresp<-max(alldata$response)
		minresp<-min(alldata$response)
		rangeresp<-maxresp-minresp
		maxob<-maxresp + rangeresp*length(unique(levels(as.factor(alldata$Type))))*0.075
		minob<-minresp - rangeresp*0.1

		cat<-c(as.factor(alldata$Type))
		catlab<-c(levels(alldata$response))

		if(is.numeric(alldata$Type)==FALSE)
		{
			par(las=1)
			plot(as.numeric(alldata$response)~as.numeric(alldata$logdose), type = "p", data=catplotdata, 	col=cat,pch=cat, xlab=XAxisTitle, ylab = YAxisTitle, ylim=c(minob,maxob))
			legend("topright", legend=levels(as.factor(tempdata$Type)),cex=0.6,pch=c(1:nlevels), bg="white", col=c(1:nlevels))
			av<-seq(-10,10,0.01)
			bv<-predict(dosefit,list(logconczzzz=av))
			lines(av,bv,col="black")
		} 
	}
	void<-HTMLInsertGraph(GraphFileName=sub("[A-Z0-9a-z,:,\\\\]*App_Data[\\\\]","", ncscatterplot3), Align="centre")

	if (DoseTransform == "Loge")
	{
		QCs$backtran<-C-(1/B)*log(((A-D)/(QCs$response-D))-1)
		length <- length(unique(levels(as.factor(QCs$dose))))
		vectorbtmean <-c(1:length)
		vectormean <-c(1:length)
		vectorStDev <-c(1:length)
		vectorN <-c(1:length)
		for (i in 1:length)
		{
			sub<-subset(QCs, QCs$dose == unique(levels(as.factor(QCs$dose)))[i])
			sub2<-data.frame(sub)
			vectorbtmean[i]=exp(mean(sub2$backtran, na.rm=TRUE))-offset
			vectorStDev[i]=sqrt((exp(2*mean(sub2$backtran, na.rm=TRUE))) * (exp(sd(sub2$backtran, na.rm=TRUE)*sd(sub2$backtran, na.rm=TRUE)))*(exp(sd(sub2$backtran, na.rm=TRUE)*sd(sub2$backtran, na.rm=TRUE))-1));

			vectormean[i]=mean(sub2$dose, na.rm=TRUE)
			tempy<-na.omit(sub2$backtran)
			vectorN[i]=length(tempy)
		}
	}

	if (DoseTransform == "Log10")
	{
		QCs$backtran<-C-(1/B)*log10(((A-D)/(QCs$response-D))-1)

		length <- length(unique(levels(as.factor(QCs$dose))))
		vectorbtmean <-c(1:length)
		vectormean <-c(1:length)
		vectorStDev <-c(1:length)
		vectorN <-c(1:length)
		for (i in 1:length)
		{
			sub<-subset(QCs, QCs$dose == unique(levels(as.factor(QCs$dose)))[i])
			sub2<-data.frame(sub)
			vectorbtmean[i]=10**(mean(sub2$backtran, na.rm=TRUE))-offset
			vectorStDev[i]=sqrt((10**(2*mean(sub2$backtran, na.rm=TRUE))) * (10**(sd(sub2$backtran, na.rm=TRUE)*sd(sub2$backtran, na.rm=TRUE)))*(10**(sd(sub2$backtran, na.rm=TRUE)*sd(sub2$backtran, na.rm=TRUE))-1));
			vectormean[i]=mean(sub2$dose, na.rm=TRUE)
			tempy<-na.omit(sub2$backtran)
			vectorN[i]=length(tempy)
		}
	}

	tempdata<-data.frame(cbind(vectormean,vectorbtmean,vectorStDev,vectorN))
	tempdata$RE<-100*(vectorbtmean-vectormean)/vectormean
	tempdata$CV<-100*vectorStDev/vectorbtmean
	colnames(tempdata)<-c("True QC mean", "Backtransformed QC mean", "Standard deviation of QC samples", "No. of backtransformable QCs","Relative error (%)", "Coefficient of variation (%)")

	HTMLbr()
	title<-c("Table of QC summary statistics")
	HTML.title(title, HR=2, align="left")

	ablen<-length(unique(rownames(tempdata)))
	tab<-matrix(nrow=ablen, ncol=11)

	for (i in 1:ablen)
	{
		tab[i,1]=format(round(tempdata[i,1],2), nsmall=2, scientific=FALSE)
	}
	for (i in 1:ablen)
	{
		tab[i,2]=c(" ")
	}
	for (i in 1:ablen)
	{
		tab[i,3]=format(round(tempdata[i,2], 2), nsmall=2, scientific=FALSE)
	}
	for (i in 1:ablen)
	{
		tab[i,4]=c(" ")
	}
	for (i in 1:ablen)
	{
		tab[i,5]=format(round(tempdata[i,3], 3), nsmall=3, scientific=FALSE)
	}
	for (i in 1:ablen)
	{
		tab[i,6]=c(" ")
	}
	for (i in 1:ablen)
	{
		tab[i,7]=format(round(tempdata[i,4], 0), nsmall=0, scientific=FALSE)
	}
	for (i in 1:ablen)
	{
		tab[i,8]=c(" ")
	}
	for (i in 1:ablen)
	{
		tab[i,9]=format(round(tempdata[i,5], 2), nsmall=2, scientific=FALSE)
	}
	for (i in 1:ablen)
	{
		tab[i,10]=c(" ")
	}
	for (i in 1:ablen)
	{
		tab[i,11]=format(round(tempdata[i,6], 2), nsmall=2, scientific=FALSE)
	}
	
	colnames(tab)<-c("True QC mean", "   ", "Back calculated QC mean", "   ","St dev of back calculated QC mean", "   ", "No. of back calculated QCs", "  ", "Relative error (%)", "   ", "Coefficient of variation (%)")

	HTML(tab, classfirstline="second", align="left")

	add4<-c("Note: The relative error (%) and coefficient of variation (%) for an individual QC are only reliable statistics if all the QCs can be back calculated. ")

	HTML.title("</bf> ", HR=2, align="left")
	HTML.title(add4, HR=0, align="left")
}


#Samples

if (Samples  != "NULL" )
{
	
	HTMLbr()
	title<-c("Back calculated sample responses")
	HTML.title(title, HR=2, align="left")

	statdata$samplesresponsezzzz = eval(parse(text = paste("statdata$", Samples)))


	if (DoseTransform == "Loge")
	{
		statdata$backtranss<-format(round(exp(C-(1/B)*log(((A-D)/(statdata$samplesresponsezzzz-D))-1))-offset, 3), nsmall=3, scientific=FALSE)
	}
	if (DoseTransform == "Log10")
	{
		statdata$backtranss<-format(round(10**(C-(1/B)*log10(((A-D)/(statdata$samplesresponsezzzz-D))-1))-offset, 3), nsmall=3, scientific=FALSE)
	}

	for (i in 1:length(statdata$samplesresponsezzzz))
	{
		statdata$blank[i]<-c(" ")	
	}


	samples<-cbind(statdata$samplesresponsezzzz, statdata$blank,statdata$backtranss)

	row<-c(" "," "," ")
	samples2<-rbind(row, samples)
	samples2<-na.omit(samples2) 

	samlen<-dim(samples2)[1]

	index<-c("Sample ID",1:(samlen-1))
	rownames(samples2)<-index
	colnames(samples2)<-c("    Sample response    ", "      ","    Back calculated response    ")
	
	HTML(samples2, classfirstline="second", align="left")
}

#end of Fourparamter if statement
}

#References
HTMLbr()
HTML.title("<bf>R references", HR=2, align="left")

HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>   R Development Core Team (2008). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.", HR=0, align="left")

#mtvnorm
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>   Alan Genz, Frank Bretz, Torsten Hothorn with contributions by Tetsuhisa Miwa, Xuefei Mi, Friedrich Leisch and Fabian Scheipl	(2008). mvtnorm: Multivariate Normal and t Distributions. R package version 0.9-0.
	", HR=0, align="left")

#lattice
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf> 
	Deepayan Sarkar (2009). lattice: Lattice Graphics. R package version 0.17-22. http://CRAN.R-project.org/package=lattice
	", HR=0, align="left")

#R2HTML
HTML.title("<bf> ", HR=2, align="left")
HTML.title("<bf>
	Lecoutre, Eric (2003). The R2HTML Package. R News, Vol 3. N. 3, Vienna, Austria.
	", HR=0, align="left")


quit()


dosefit<-nls(bweight2~D + (A-D)/(1+10^((C-logdosenumber)*B)), start=list(A=maxp,B=slopep,C= ed50p,D=minp), data=statdata)

par(mfrow = c(1,2))
plot(logdosenumber,statdata$bweight2, main = "Non linear fit�)
av<-seq(-10,10,0.01)
bv<-predict(dosefit,list(logdosenumber=av))
lines(av,bv,col=�red�)
summary(dosefit)
plot(fitted(dosefit), residuals(dosefit),main="Pred vs. Resid")